#include "libRNG.h"
#include "oRNGs.h"

#pragma comment(lib,"cudart.lib")
#pragma comment(lib,"curand.lib")

#pragma comment(lib,"mkl_core.lib")
#pragma comment(lib,"mkl_intel_c.lib")
#pragma comment(lib, "mkl_intel_thread.lib")

void createRNG( RNG_ID method, RNGStatus *status, unsigned int seed )
{
	status->method = method;

	switch (method)
	{
	case CURAND1:		
		curandCreateGenerator( &status->gen, CURAND_RNG_PSEUDO_XORWOW );		
		curandSetPseudoRandomGeneratorSeed( status->gen, (unsigned long long)seed );			
		break;
	case CURAND2:		
		curandCreateGenerator( &status->gen, CURAND_RNG_PSEUDO_XORWOW );		
		curandSetPseudoRandomGeneratorSeed( status->gen, (unsigned long long)seed );		
		curandSetGeneratorOrdering ( status->gen, CURAND_ORDERING_PSEUDO_SEEDED );
		break;
	case CURAND3:		
		curandCreateGenerator( &status->gen, CURAND_RNG_PSEUDO_MRG32K3A );		
		curandSetPseudoRandomGeneratorSeed( status->gen, (unsigned long long)seed );			
		break;
	case CURAND4:		
		curandCreateGenerator( &status->gen, CURAND_RNG_PSEUDO_MTGP32 );		
		curandSetPseudoRandomGeneratorSeed( status->gen, (unsigned long long)seed );			
		break;
	case CURAND5:		
		curandCreateGenerator( &status->gen, CURAND_RNG_PSEUDO_PHILOX4_32_10 );		
		curandSetPseudoRandomGeneratorSeed( status->gen, (unsigned long long)seed );			
		break;		

	case MKL1: 
		vslNewStream( &status->stream, VSL_BRNG_MCG31, seed ); 
		break;
	case MKL2:
		vslNewStream( &status->stream, VSL_BRNG_R250, seed ); 
		break;
	case MKL3:
		vslNewStream( &status->stream, VSL_BRNG_MRG32K3A, seed ); 
		break;
	case MKL4:
		vslNewStream( &status->stream, VSL_BRNG_MCG59, seed ); 
		break;
	case MKL5:
		vslNewStream( &status->stream, VSL_BRNG_WH, seed ); 
		break;
	case MKL6:
		vslNewStream( &status->stream, VSL_BRNG_MT19937, seed ); 
		break;	
	case MKL7:
		vslNewStream( &status->stream, VSL_BRNG_MT2203, seed ); 
		break;	
	case MKL8:
		vslNewStream( &status->stream, VSL_BRNG_SFMT19937, seed ); 
		break;		

	case LCG:
		srand(seed); 		
		vslsNewAbstractStream( &status->stream, BUFN, sbuf, 0, 1, sLCG );	
		vsldNewAbstractStream( &status->stream_d, BUFN, dbuf, 0, 1, dLCG );	
		break;	

	case KISS:
		kiss_srand((unsigned long)seed); 		
		vslsNewAbstractStream( &status->stream, BUFN, sbuf, 0, 1, sKISS );
		vsldNewAbstractStream( &status->stream_d, BUFN, dbuf, 0, 1, dKISS );
		break;
	}
}
void destroyRNG( RNGStatus *status )
{
	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		curandDestroyGenerator( status->gen );			
		break;

	case LCG:
	case KISS:
		vslDeleteStream( &status->stream );		
		vslDeleteStream( &status->stream_d );
		break;
		
	default:
		vslDeleteStream( &status->stream );		
		break;	
	}
	status->method = ILLEGAL;
}

template<typename T>
__global__ void scale_kernel( T * vector, T alpha, T beta, int count )
{
	for ( int tidx =  threadIdx.x + blockIdx.x*blockDim.x;
		tidx < count;
		tidx +=  blockDim.x*gridDim.x )
	{
		vector[tidx] = alpha*vector[tidx]+beta;
	}
}

template<typename T>
void scale( T* vector, T alpha, T beta, int count )
{
	int nthreads = 32*2;
	int nblocks = 8*3;
	scale_kernel<<<nblocks,nthreads>>>( vector, alpha, beta, count );
}

void generateUniform( RNGStatus *status, float* randoms, int count, float lower, float upper )
{	
	float* devData;
	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		cudaMalloc( &devData, sizeof(float)*count );
		curandGenerateUniform( status->gen, devData, count );
		scale( devData, upper-lower, lower, count );
		cudaMemcpy( randoms, devData, count * sizeof(float), cudaMemcpyDeviceToHost );			
		cudaFree( devData );
		break;

	default:
		vsRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream, count, randoms, lower, upper );			
		break;
	}
}
void generateUniformDouble( RNGStatus *status, double* randoms, int count, double lower, double upper )
{	
	double* devData;
	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		cudaMalloc( &devData, sizeof(double)*count );
		curandGenerateUniformDouble( status->gen, devData, count );
		scale( devData, upper-lower, lower, count );
		cudaMemcpy( randoms, devData, count * sizeof(double), cudaMemcpyDeviceToHost );	
		cudaFree( devData );
		break;

	case LCG:
	case KISS:
		vdRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream_d, count, randoms, lower, upper );
		break;

	default:
		vdRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream, count, randoms, lower, upper );			
		break;
	}
}

void generateNormal( RNGStatus *status, float* randoms, int count, float mean, float stdev )
{	
	float* devData;
	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		cudaMalloc( &devData, sizeof(float)*count );

		curandGenerateNormal( status->gen, devData, count, mean, stdev );
		cudaMemcpy( randoms, devData, count * sizeof(float), cudaMemcpyDeviceToHost );

		cudaFree( devData );
		break;

	default:
		vsRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream, count, randoms, mean, stdev );		
		break;	
	}	
}
void generateNormalDouble( RNGStatus *status, double* randoms, int count, double mean, double stdev )
{	
	double* devData;
	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		cudaMalloc( &devData, sizeof(double)*count );

		curandGenerateNormalDouble( status->gen, devData, count, mean, stdev );
		cudaMemcpy( randoms, devData, count * sizeof(double), cudaMemcpyDeviceToHost );

		cudaFree( devData );
		break;

	case LCG:
	case KISS:
		vdRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream_d, count, randoms, mean, stdev );
		break;
	default:
		vdRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream, count, randoms, mean, stdev );		
		break;
	}	
}

void generateUniform_d( RNGStatus *status, float* randoms, int count, float lower, float upper )
{	
	float* hostData;	

	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		curandGenerateUniform( status->gen, randoms, count );
		scale( randoms, upper-lower, lower, count );
		break;

	default:
		hostData = (float*)malloc( sizeof(float)*count );
		vsRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream, count, hostData, lower, upper );	
		cudaMemcpy( randoms, hostData, sizeof(float)*count, cudaMemcpyHostToDevice );
		free( hostData );
		break;
	}		
}
void generateUniformDouble_d( RNGStatus *status, double* randoms, int count, double lower, double upper )
{	
	double* hostData;	

	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		curandGenerateUniformDouble( status->gen, randoms, count );
		scale( randoms, upper-lower, lower, count );
		break;

	case LCG:
	case KISS:
		hostData = (double*)malloc( sizeof(double)*count );
		vdRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream_d, count, hostData, lower, upper );	
		cudaMemcpy( randoms, hostData, sizeof(double)*count, cudaMemcpyHostToDevice );
		free( hostData );
		break;	


	default:
		hostData = (double*)malloc( sizeof(double)*count );
		vdRngUniform( VSL_RNG_METHOD_UNIFORM_STD, status->stream, count, hostData, lower, upper );	
		cudaMemcpy( randoms, hostData, sizeof(double)*count, cudaMemcpyHostToDevice );
		free( hostData );
		break;	
	}

}

void generateNormal_d( RNGStatus *status, float* randoms, int count, float mean, float stdev )
{	
	float* hostData;

	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		curandGenerateNormal( status->gen, randoms, count, mean, stdev );		
		break;

	default:
		hostData = (float*)malloc( sizeof(float)*count );
		vsRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream, count, hostData, mean, stdev );
		cudaMemcpy( randoms, hostData, sizeof(float)*count, cudaMemcpyHostToDevice );
		free(hostData);
		break;	
	}		
}
void generateNormalDouble_d( RNGStatus *status, double* randoms, int count, double mean, double stdev )
{	
	double* hostData;

	switch ( status->method )
	{
	case CURAND1:			
	case CURAND2:			
	case CURAND3:			
	case CURAND4:
	case CURAND5:
		curandGenerateNormalDouble( status->gen, randoms, count, mean, stdev );		
		break;

	case LCG:
	case KISS:
		hostData = (double*)malloc( sizeof(double)*count );
		vdRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream_d, count, hostData, mean, stdev );
		cudaMemcpy( randoms, hostData, sizeof(double)*count, cudaMemcpyHostToDevice );
		free(hostData);
		break;	

	default:
		hostData = (double*)malloc( sizeof(double)*count );
		vdRngGaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER2, status->stream, count, hostData, mean, stdev );
		cudaMemcpy( randoms, hostData, sizeof(double)*count, cudaMemcpyHostToDevice );
		free(hostData);
		break;	
	}		
}
